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Abstract 



To understand the origin of ultra-high energy cosmic rays (UHECRs, defined to 
be above 10^^ eV), it is required to model in a realistic way their propagation in 
the Universe. UHECRs can interact with low energy radio, microwave, infrared and 
optical photons to produce electron/positron pairs or pions. The latter decay and 
give rise to neutrinos and electromagnetic cascades extending down to MeV ener- 
gies. In addition, deflections in cosmic magnetic fields can influence the spectrum 
and sky distribution of primary cosmic rays and, due to the increased propagation 
path length, the secondary neutrino and 7— ray fluxes. Neutrino, 7— ray, cosmic ray 
physics and extra-galactic magnetic fields are, therefore, strongly linked subjects 
and should be considered together in order to extract maximal information from 
existing and future data, like the one expected from the Auger Observatory. For 
that purpose, we have developed CRPropa, a publicly-available numerical package 
which takes into account interactions and deflections of primary UHECRs as well as 
propagation of secondary electromagnetic cascades and neutrinos. CRPropa allows 
to compute the observable properties of UHECRs and their secondaries in a vari- 
ety of models for the sources and propagation of these particles. Here we present 
physical processes taken into account as well as benchmark examples; a detailed 
documentation of the code can be found on our web site. 
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1 Introduction 



Astroparticle physics is currently experimentally driven and involves many dif- 
ferent existing or planned projects ranging from UHECR observatories such as 
the Pierre Auger Observatory [1], to neutrino telescopes [2], as well as ground 
and space based 7— ray detectors operating at TeV and GeV energies, respec- 
tively [3]. It is clear that GeV- TeV 7— ray and neutrino astronomy will prove 
an invaluable tool to unveil the sources, and probe into the mechanism, of 
UHECRs. Even if a putative source were to produce exclusively UHECRs, 
photo-pion [4] and pair production by protons on the cosmic microwave back- 
ground (CMB) would lead to guaranteed secondary photon and neutrino fluxes 
that could be detectable. Furthermore, spectra, power and sky distributions 
of both primary UHECRs and secondary 7— rays and neutrinos depend on the 
poorly known large scale cosmic magnetic fields. 

It is, therefore, desirable to have a numerical tool that can treat the interface 
between UHECR, 7— ray and neutrino astrophysics, and large scale magnetic 
fields. To this end, we have recently merged our Monte Carlo code for simu- 
lating three dimensional propagation of UHECRs in a structured, magnetized 
Universe [5] with a one- dimensional transport code that solves electromagnetic 
(EM) cascades and neutrino propagation [6]. We discuss the limitations due 
to the one-dimensional approximation and implement a procedure to test the 
resulting uncertainty on the EM cascade on the observable fluxes. With the 
present paper, we release a public version of this code which we hope will be 
useful for the cosmic ray, 7— ray and neutrino communities. 

In the following, we present the relevant interactions and propagation phenom- 
ena taken into account, and the propagation algorithms applied in CRPropa. 
We also present a few examples of how to use the code in practice. The nu- 
merical package and its detailed documentation are available for downloading 
on the CRPropa website, ,http : / / apcauger . in2p3 . f r/ CRPropa, . 

We use natural units, c = h = 1 throughout this paper. 



2 Propagation algorithms 

UHECRs are injected at specified sources, and propagated step-by-step in ei- 
ther a one- or a three-dimensional environment. The trajectories are regularly 
sampled, or recorded only at specific locations (e.g. at a given distance from 
a source, or at an "observer" point). Each propagation step consists of inte- 
grating the Lorentz equations, and computing the interactions and possibly 
the secondaries generated by those interactions. 
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Fig. 1. Principle of the propagation algorithm. This scheme applies to all configu- 



In the 3-dimensionnal "simulation box" is defined and periodic bound- 

ary conditions are assumed. 

When deflections are taken into account, cosmological redshifts cannot be 
computed, because the propagation time until the particle reaches the observer 
is not known before hand. Therefore, redshift evolution is only accounted for in 
the ID version of the package. The concordance cosmology is used for which, 
assuming a flat Universe, the Hubble rate ii{z) at redshift z in the matter 
dominated regime, z < 10^, is given by 



The parameters and Oa can be freely chosen, their standard values being 
= 0.3, riA = 0.7, and Hq = ho 100 km s"^ Mpc"^ with ho = 0.72. 

The general principle of the simulations is shown in Fig. 1. 

2.1 Nucleon Interactions 

The most famous interaction of nucleons with the low-energy photon back- 
grounds is pion production, which generates the GZK feature. In order to 
handle pion production, we use the event generator SOPHIA [7], that has 
been explicitely designed to study this phenomenon and that uses the particle 
production cross-sections measured in accelerators. We have also augmented 
the SOPHIA package for interactions with a low energy extragalactic back- 
ground light (EBL) with a general energy distribution. SOPHIA allows to 
determine the distribution of the stables particles generated by an interaction 
with a low-energy photon. 
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Pair production by protons (PPP) on the CMB, also known as Bethe-Heitler 
process, is taken into account as a continuous energy loss whose rate we eval- 
uate following the expressions in Refs. [8,9]. For the spectrum of the pairs we 
exploit the fact that Bethe-Heitler and triplet pair production, e^b ~^ ee'^e~ , 
are analogous electromagnetic processes, their cross sections and inelasticities 
converging for relativistic pairs. Fig. 2 of Ref. [10] then shows that the spec- 
trum of electron-positron pairs (heretoafter simply referred to as electrons) 
generated by a proton of energy E can be approximated by a power-law en- 
ergy distribution dn/dEg oc E'^'^/^. Kinematics implies that this power law 
holds for E'min ^ Eg < Eppp, where the minimal and maximal energies are 
given by [6] 



E ^ 4.5xlO^^(^)^(^)eV 
+ 4.6 X 10- (^) (^) -h 1 



__3^_3.3xlO-^^J eV. (2) 

In Eq. (2), rrip and rrie are the proton and electron masses, respectively, 
e is the low energy target photon energy, and the approximation for -Emm 
holds for rriemp < AEe < m^. The average electron energy is then E^. = 

/Ij/ dEgEgEf/^l dEeE-'^'^ ~ 3 E^^^^E^pp which is indeed much smaller 
than the primary proton energy E. Prom Eq. (2), the inelasticity K = E^/E, 
whose precise energy dependence can be found in Ref. [9] , for meirip < AEe < 
mi can thus be approximated by 



3 m^/^ 



~3.4 X 10"^ 



/ E e ^-V2 

VEeVy 



meV 



This is consistent with Figs. 1 and 2 in Ref. [11]. For our purposes, we are 
not sensitive to the lower kinematic limit since the total energy produced oc 
fi^^f dEgEgEg'^/'^ ~ 4£'ppp is insensitive to E'min as long as E'min ^ -Eppp, but 
rather is dominated by the highest energies. As a consequence, the total proton 
energy loss rate due to pair production is dominated by the highest energy 
electrons close to Eppp. However, because the production cross section of these 
highest energy electrons is much smaller than the one for the more numerous 
lower energy electrons, the average inelasticity Eq. (3) is nevertheless small, 
below 10~^ everywhere above the pair production threshold. The spectrum and 
maximal energy of the pairs will be important for the synchrotron spectrum 
emitted by these electrons in an EGMF of strength B which peaks at ~ 
6.8 X lO^i (Ee/10i9eV)2(E/0.1/iG)eV. 
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Nuclcons can be followed down to 10^ eV with CRPropa, below which inter- 
actions become negligible. 



2.2 Secondary Electromagnetic Cascades and Neutrinos 



The secondary neutrinos from pion production of nucleons are propagated in 
straight lines assuming no energy losses except redshift effects. 

All the EM products of these interactions are evolved using an EM cascade 
code based on Ref. [6]. The photons and pairs are followed until either their en- 
ergy drops below 100 MeV or they reach an observer. All relevant interactions 
with a background photon 7b are taken into account, namely single pair pro- 
duction (PP), 776 e+e^, double pair production (DPP), 775 — > e+e^e+e", 
inverse Compton scattering (ICS), e7b 67, and triplet pair production 
(TPP), e7b — > ee+e~ (see also Ref. [12] for a detailed discussion of imple- 
mented interactions). In addition, synchrotron losses of electrons in the (in 
general) inhomogeneous EGMF are taken into account and the resulting lower 
energy synchrotron photons are also followed in the subsequent EM cascade. 

This module has been applied to EM cascades from discrete magnetized proton 
sources in galaxy clusters in Ref. [13] . The EM cascades that are followed with 
the current version of CRPropa are propagated in straight lines, even in the 
case of 3-dimensionnal simulations for UHECRs: Every time a primary hadron 
interacts and initiates an EM cascade, it is assumed that the secondaries 
propagate along straight lines and it is checked whether the line of sight crosses 
the observer. If this is the case, the EM cascade module is called with the 
corresponding propagation distance and the projected magnetic field profile. 
Electrons in the EM cascade can of course be deflected in the EGMF, and we 
discuss here the validity of this one-dimensionnal approximation. 

In a magnetic field of strength B the synchrotron cooling time for an electron 
of energy is given by 



E. \-WB\-' 



~3.84kpc 



1015 eV; \iiG^ 



where ctt = 87rQ;^/3mg is the Thomson cross section, with a the fine structure 
constant. At high energies, in the Klein-Nishina regime the inverse Compton 
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energy loss length is roughly [12] 



iic < 400pc 



) 



for£;e > lO^'eV. 



(5) 



1015 eV 



At energies Eg > 10^^ eV in Eq. (5) the energy loss length is between a factor 
~ 30 and a few hundred smaller than the numerical value in Eq. (5) due 
to contributions from the universal radio background. For a conservatively 
large tic at these energies wc use an interpolation of Fig. 12 in [12] for the 
conservatively low radio background estimate. For < 10^^ eV, ICS on the 
CMB is in the Thomson regime, with an interaction length Ajc ~ 1/o"tT^cmb ~ 
1.2 kpc, with ncMB the CMB photon density. The energy lost by an electron at 
each interaction is 5£'e ~ AeE^/Sml, where e is a typical CMB photon energy. 
As a consequence, the energy loss length at energies below ~ 10^^ eV is: 



These length scales as well as the maximal propagation distance must be 
compared with the Larmor radius 



In order for a one-dimensional treatment of EM cascades to be a good ap- 
proximation, the Larmor radius has to be much larger than either the to- 
tal propagation length, the IC or the synchrotron loss lengths. For a given 
magnetic field, the condition ri > A x min(tsynch, tic) results in a condition 
Ee > Ec{B,A), corresponding to deflections of the pairs by < {10/ A) x 6°. 
This estimate of the deflection angle is however conservatively large: in real- 
istic situations, magnetic fields are inhomogeneous with many reversals along 
the fine of sight, and the actual defiection angle will be smaller provided the 
magnetic field coherence length is smaller than the energy loss length. The 
dependence of the "critical" energy Ec on B and A is shown in Fig. 2: For 
yl = 10, corresponding to defiection by less than ^ 6°, Ec is determined by the 
competition between defiections and ICS for B < 300 pG, or between defiec- 
tion and synchrotron emission for B > 300 pG. For A — 100, corresponding 
to defiection by less than ~ 0.6°, the transition between ICS and synchrotron 
emission as dominant losses to be compared with defiection occurs at S ~ 20 



It turns out that whenever Ec{B, A) < 10^^ eV, the 7— ray fiux from defiected 
pairs is sub-dominant. This is simply due to the fact that the pair fiux at ener- 
gies Eg < 10^5 eV is suppressed compared to the 7— ray fiux which has a much 
larger interaction length and piles up below the pair production threshold. The 
7— ray fiux from defiected pairs can only be important if Ec{B, A) ^ 10^^ eV 







pG. 
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Fig. 2. The critical energy E^, below which the e^/~ are deflected before cascading 
to lower energies, as a function of the order of magnitude of the magnetic field. Ec 
is obtained with the parameterization of various timescales given in the text, for 
A = 10 (solid line) and A = 100 (dashed line). This corresponds to cutting all pairs 
being deflected by more than ~ 6° or 0.6°, respectively. Note that the jump around 
~ 3 X lO"^'^ G and ~ 2 x 10^^^ G, respectively, is due to the transition from ICS to 
synchrotron emission (at large fields) as dominant energy loss. 

which, from Fig. 2, requires that synchrotron emission dominates the losses of 
the deflected pairs. In this case, a significant fraction of the energy flux going 
into pairs is defiected more than ~ (10/^) x 6°, thus modifying the 7— ray 
point fiux at energies 

.,...x.O«ev(^)^(A). 

In the foUowing we wiU confirm these expectations with numerical simulations. 

Within CRPropa, the parameter A can be chosen by the user, and the local 
contribution of electrons with energy Eg < Ec{B,A) to the 7-ray fiux can be 
switched on or off: This allows to estimate the uncertainty in the 7-ray flux 
arriving within a certain angle ~ (lO/A) x 6° from a point source due to the 
ID approximation. An example is shown in Fig. 3, where the computed 7-ray 
fluxes from a single proton source located at 100 Mpc from the Earth in a 
uniform magnetic fleld of amplitude 100 pG are compared with and without 
cutting the charged component of the EM cascade deflected by more than 6° 
and 0.6°, respectively. 

For the flux arriving within 6° in Fig. 3, Ec — S x 10^^ eV, see Fig. 2, and 



7 




Fig. 3. Flux of secondary neutrinos (all flavors are added) and 7-rays from a single 
source of UHE protons with injection spectrum oc E~'^ up to 5 x 10^'^ eV, computed 
assuming a straight line propagation for the protons but taking into account the 
influence of a 100 pG magnetic field on the EM cascades. The red line is the flux 
computed by "cutting" the local e+/~ flux below Ec{B,A) (see text), for A = W 
(continuous line) and A = 100 (dashed line). This corresponds to pair deflections of 
6° or 0.6°, respectively. 

indeed a discernible but still modest, ~ 30%, modification appears only for 
Ey < 0.1 TeV, where the photon energy flux becomes comparable to the pair 
energy flux around Ec- 

For the flux arriving within 0.6° in Fig. 3, E^ — 2 x 10^^ eV, see Fig. 2, and by 
the above argument and Eq. (8) we expect the photon flux to be significantly 
modified below ~ 10 TeV. Indeed, at these energies the flux is reduced by a 
factor ~ 5. 

In case of comparatively strong magnetic fields of order fiG, typical in galaxy 
clusters, Ec < 10^^ eV, and 7— ray point fluxes arriving within ~ 0.6° should 
only be modified significantly for E^ < 100 GeV. Also note that for the pro- 
duction of secondaries inside a magnetized region where the parent UHECR 
particles are isotropically distributed, the full three dimensional treatment of 
the EM cascade is not necessary because for any e~ that is deflected away 
from the line of sight there is always another e~ that is deflected into the 
line of sight. In realistic situations, the magnetic fields are highly structured 
with typical amplitudes of ~ fiG in the clusters, and < 10 pG in the voids. 
The above discussion shows that in all these cases CRPropa can estimate the 
minimal 7— ray flux arriving within an angle (lO/A) x 6° from the source. 
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Fig. 4. Models implemented for the low energy photon background at zero redshift. 
The IRB consists basically of a peak in the far infrared around 100/im dominated 
by dust and a peak in the near infrared dominated by stars. 

2.3 Background Photon Spectra and their Evolution 

Fig. 4 shows the EBL energy distributions that have been implemented. The 
most important is the CMB. For the infrared background (IRB) we imple- 
mented three distributions, a low and a high version of Franceschini et al. [14] 
which differ roughly by a factor 5, as well as the one by Primack et al. [15]. The 
low Franceschini et al. and the Primack et al. backgrounds are consistent with 
recent upper limits from blazar observations in TeV 7— rays by HESS [16]. For 
a recent review of the IRB see for example Ref. [17]. 

The IRB has a significant influence on EM cascades only around the threshold 
for pair production, i.e. between a few Tev and ~ 100 TeV. At higher ener- 
gies, the 7— ray flux is suppressed by interactions with the CMB and, above 
~ 10^^ eV, by interactions with with the radio background. At energies below 
~ TeV, the Universe acts as a calorimeter and the total photon flux is propor- 
tional to the total EM energy injected above ~ PeV with a rather universal 
shape [18]. 

Although its photon number density ^ 2 cm~^ is a factor ~ 200 smaller than 
for the CMB, below the GZK-cutoff and above ~ 10^'' eV the IRB can signif- 
icantly reduce the nucleon mean free path for pion production. This can be 
important for secondary photon and neutrino [19,20] production, especially 
for a steep primary injection spectrum and/or strong redshift evolution. 
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Fig. 5. Proton energy loss length for pair production on the CMB (continuous line), 
interaction length for pion production on the CMB (dashed line) and on the Primack 
et al. IRE (dotted line) at 2; = 0. The irregularities in the dashed curve are due to 
the piecewise power law fits of the Primack et al. IRB. 

For the universal radio background (URB) we use a weak and a strong version 
based on Rcf. [21] and on observations [22]. The URB is mostly important for 
EM cascades above ~ 10^^ eV where it can inhibit cascade development due 
to the resulting small pair production lengths, especially for fast synchrotron 
losses of electrons in the presence of strong magnetic fields. 

Since URB photons can give rise to pion production only above a few times 
10^^ eV, where the interaction rate is essentially proportional to the total EBL 
photon density which is dominated by the CMB by a factor ~ 10^, see Fig. 4, 
the URB is negligible for pion production. The same applies to pair production 
by protons. 

Figs. 5 and 6 show interaction and energy loss lengths for protons and inter- 
action lengths of photons, respectively, and their dependence on EBL models 
at zero redshift. This demonstrates that the IRB becomes important for pion 
production by protons below the GZK cutoff and for pair production by pho- 
tons below the threshold in the CMB at ~ 10^^ eV. It also shows that the 
URB tends to dominate pair production by photons above ~ 10^^ eV. 

The redshift evolution of the cosmic microwave background (CMB) is trivial. 

The redshift evolution of the radio and infrared distributions is more compli- 
cated: Ultra-relativistic particles of energy E injected at redshift z' with a rate 
per energy and comoving volume $(£^, z') result in a physical number density 
per energy at redshift z given by 
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Fig. 6. Photon interaction length at z = on the EBL consisting of the CMB, the 
Primack et al. IRB, and the strong URB version. Dotted hne: Interaction length in 
the CMB only at z = 0. 



n{E,z)^{l + zf [ 

J z 



dz' 



,47r$ [E,{E,z,z'),z'] 



{l + z')H{z') 

x^{E,z,z'), (9) 

where it is assumed that the particle looses energy continuously such that its 
injection energy can be computed analytically, Ei{E, z, z'). Interactions of the 
low energy EBL photons, whose differential number densities we will denote 
by nf,(e, z) in the following to distinguish from the high energy particles, can 
safely be neglected after recombination, z < 10^, such that Ei{E,z,z') = 
(1 + z')E/{l + z). Eq. (9) then simplifies to 

= + dJ ^ (10) 

By using \dtldz\ — [(1 + z)B.{z)\~^ ^ one can see easily that the total en- 
ergy density per comoving volume redshifts as / deenb{s, z)/{l + zY = (1 + 
z) J dt dsi ^{Si, z') / {1 + z') , as it should be. 

For the URB we implemented a nontrivial redshift evolution in the cascade 
module, as this can be relevant for EM cascade development. We assume that 
^vkb{.^iZ) = 0urb(£)5'urb(-2) factorizes into an energy dependence 0urb(£) 
motivated by the observations [22] and theoretical estimates [21] and a redshift 
dependence given by 

^urb(^) = 10^-^«-°-^«^^ , (11) 
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as in Ref. [6]. 



For the Primack et al. IRB [15] we use for simplicity the differential photon 
energy distribution evolution 



nb{e,z) 



(1 + zfrib (j^, ^ = O) ior z<Zb, 
otherwise 



;i2) 



which corresponds to instantaneous creation of the background at redshift 
Zb with = H{zb)nb[e/{1 + Zb),z = 0]6{z' - Zb)/{4n) in Eq. (10). It 

strictly applies to the CMB which was effectively produced at decoupling, 
Zb ~ 1100. For the IRB we assume Zb — 5. Interaction lengths l{E,z) and, 
in case of continuous energy loss processes such as PPP, energy loss rates 
h{E, z) = dE/dt then follow simple scaling relations in redshift [20], 



1{E, z)-^ = (1 + zfl [(1 + z)E, z = 0]"' 
h{E, z) = (1 + z)% [(1 + z)E, z = 0] . (13) 

This simplifies implementation in SOPHIA. 



2.4 Distributions and Properties of Sources 



Both single sources and reahzations of both discrete or continuous source 
distributions can be used in CRPropa. In the latter case, the distributions 
can be selected, for example, to follow the baryon density from a large scale 
structure simulation box, and are periodically repeated. 

The UHECR particles are injected isotropically around the sources with a 
monochromatic or a power-law energy distribution between a minimal and a 
maximal energy, E'min and fi'max, respectively: 



oc E ■ < E- ■ < E 

For each trajectory reaching the observer and being registered, the source 
identity i is also registered. This allows to apply a re-weighting procedure on 
the recorded "events", in order to vary individual source properties such as 
their injection power law index or luminosity Qj. For example, it is most 
efficient in terms of CPU time to inject the UHECRs with a spectral index 
cto = 1 at the sources, that is with a uniform distribution in the logarithm of 
the energy. By re- weighting each recorded event by a factor w oc QiE^^^, the 
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source i would contribute with a power Qi and an effective injection power law 
index ccj in all observables constructed from the weighted trajectory sample. 



3 Large Scale Structure and Magnetic Fields 

The strength and distribution of the EGMF is currently poorly known and 
their impact on UHECR are hard to quantify, as demonstrated by the different 
results in Refs. [5,23]. See also Ref. [24] for a discussion of these differences and 
Ref. [25] for a review on EGMF. We note that there are recent observational 
hints of EGMF as strong as ~ 0.1 fiG on scales as extended as superclus- 
ters [26], as well as theoretical motivations for such fields [27]. 

Enhanced magnetic fields around large scale structures such as galaxy clusters 
together with associated larger EBL densities can lead to increased production 
of 7— rays and neutrinos. 

The EGMF from the large scale structure simulation from Ref. [28,29] has so 
far been implemented in CRPropa, but any magnetic field model can be used. 
Within the public package CRPropa, only a small subgrid of the simulations 
from [28,29] is provided in order to allow simple tests. Fig. 7 shows a 2D cross 
section through the environment of a galaxy cluster from this simulation. In 
this simulation, the magnetic fields follow the baryon density, and in particular 
the regions that are filled with sub-//G fields are quite extended around the 
large-scale structures (with a typical extension of a few Mpc). This is due, in 
particular, to the fact that magnetic fields are generated at the LSS shocks 
within that model. Of course, the properties of 7-ray sources associated with 
UHECR sources as well as the feasibility of "charged particle astronomy" 
depend strongly on the magnetic field model [24] , [23] . 

Large scale structure simulations usually cover only a small fraction of today's 
Universe, typically of order 100 Mpc in linear scale. Since sources at much 
larger, cosmological distances can contribute to the fiuxes of UHECR below the 
GZK cutoff, of photons below ~ TeV and of neutrinos, the EGMF and source 
distributions are periodically continued in the 3D version of the code. EGMF 
with homogeneous statistical properties and power law spectra in Fourier space 
(e.g a Kolmogorov spectrum) have also been implemented in the package. 



4 Simple Applications 

We present here applications of CRPropa that are obtained with very simple 
configurations requiring little CPU time. The results can easily be compared 
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Fig. 7. A 2D cross section through the relative size and polarization of the EGMF 
in linear scaling, (top panel) and the relative baryon density in logarithmic scaling 
(bottom panel) in the environment of a galaxy cluster from the simulations from 
Ref. [28,29]. 

with previous results from the literature. 

Fig. 8 shows the averages and dispersions of the energy of nucleons in a one- 
dimensional simulation, as a function of propagated distance for various initial 
energies. Using SOPHIA automatically enables us to reproduce the stochas- 
ticity of pion production. 

Fig. 9 shows the spectra of secondaries generated during the one-dimensional 
propagation of UHECRs from a source located at 20 Mpc or 100 Mpc from 
the observer. Note that the neutrino flux increases with distance to the source, 
whereas the photon flux above ~ 10^"^ eV decreases, but the photon flux below 
this energy increases. This is because more secondary neutrinos and EM par- 
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Fig. 8. Evolution of the energy of nucleons as a function of propagation distance, for 
initial energies of 5, 50 or 500 EeV. The thin lines indicate the dispersion induced 
by the stochasticity of pion production. 

tides are produced for larger propagation distances, but EM particles above 
~ 10^^ eV are quickly degraded and cascade down to sub-PeV energies. A more 
detailed analysis of the fluxes of secondaries from a single UHECR source (e.g. 
the relative contribution of pair production and pion production on the 7-ray 
flux) can be found in Ref. [13]. The study of secondary photons from UHECR 
sources has also been carried out in various situations in Refs [30,31,32,33,34]. 

Fig. 10 shows the spectra of secondary neutrinos from a source located at 20 
Mpc from an observer, depending in particular on the magnetic field effects. 
It is remarkable that, for a given source luminosity, the flux of secondary 
neutrinos is increased by a factor of more than two due to the enhancement of 
the UHECR propagation distance generated by the /zG-level magnetic fields 
that surround this source. 

Fig. 11 compares the spectral shape of UHECRs from a source located at 100 
Mpc from an observer, depending on the presence of magnetic fields around 
the source. If magnetic fields of amplitude ~ fiG surround the source over a 
few Mpc, the observed spectrum is clearly modified: 1) there is a dispersion 
in the true propagation distance, compared to a fixed propagation distance of 
100 Mpc. This reduces the amplitude of the "bump"; 2) the mean propagation 
distance is increased compared to 100 Mpc. This leads to a GZK cut-off at 
slightly lower energies. 

Fig. 12 compares the spectra obtained with CRPropa to the one presented 
on the red curve of Fig. 14 in [35] for a model of cosmologically distributed 
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Fig. 9. Spectrum of secondary photons and neutrinos (all flavors added) generated 
by pion and pair production from a single UHECR source at a given distance. We 
consider here a one-dimensionnal model, with an injection spectral index a = 2 
for the UHECRs. A uniform magnetic field of 0.1 nG is assumed. Note that below 
'^TeV the 7— ray flux would be spread over several degrees and that, as shown in 
Fig. 3, the ID approximation of the EM cascade does not significantly affect the 
accuracy of the 7— ray flux arriving within such angles. The fluctuations at the 
highest energies are statistical. 



proton sources with spectral index 7^ = 2.6 and a source evolution parameter 
m = 2.4. We see that, for a given model, the spectrum estimations obtained 
with our Monte-Carlo method and with a direct integration of the transport 
equations (for [35]) agree within a few %. The blue and red curves of the 
lower panel in Fig. 12 show the influence of two numerical parameters on 
the accuracy of the derived spectrum at the highest energies. The maximum 
injection energy allowed in the Monte-Carlo has an influence at energies above 
10^° eV, in agreement with results shown, for example, in Fig. 5 of [35]. The 
use of a propagation stepsize of 0.3 Mpc instead of 1 Mpc does not lead 
to a significant change in the simulated spectrum. Other tests showed that 
using a propagation stepsize of 5 Mpc instead of IMpc results in a ~ 10% 
overestimation of the spectrum in the specific energy range 100 < E < 160 
EeV, and an underestimation of the spectrum at higher energies. A 1 Mpc 
stepsize is therefore reasonable to reach the typical accuracies required for 
comparison with current and forecoming experimental data. 
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Fig. 10. Secondary neutrinos (all flavors added) from a nearby source of UHECRs 
with a given luminosity. The flux increases at high energies both with maximum 
UHECR acceleration energy and with the strength of magnetic fields surrounding 
the source. The fluctuations at low energy are statistical. The y-axis is in arbitrary 
units. 




Energy (EeV) 



Fig. 11. UHECR spectrum from a source located at 100 Mpc from an observer, 
injecting protons with a spectrum oc E~'^ up to ii^max = 10^^ eV. The red curve is 
obtained from a full 3-dimensional simulation, where the source is embedded in a 
region with //G fields over a few Mpc. 
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Fig. 12. Top: Comparison of the spectra obtained with CRPropa to the one found 
in [35] for a model of cosmologically distributed sources. The specific parameters 
of the model correspond to the red curve of Fig. 14 in [35]. We use a propagation 
stepsize of IMpc and a maximum injection energy of 10^^ eV. Bottom: Relative 
difference with respect to the curve of [35] (black). The error bars are the statistical 
uncertainties due to the finite number of propagated nucleons. Red: same, for a 
simulation using a stepsize of 0.3 Mpc. Blue: same, for a simulation using a maximal 
injection energy of 10^^ eV. 

5 Conclusions 

We have presented the first pubhc package to study systematically the prop- 
erties of the propagation of UHECRs and their secondaries in a structured 
magnetized Universe. We have detailed the interactions that are already im- 
plemented, and presented a few simple examples obtained directly by running 
the CRPropa code. 

A major advantage of CRPropa is its large modularity, which should allow 
various users to implement their own modules, adapted to specific UHECR 
propagation models. Many possible upgrades of the CRPropa package can be 
considered: This includes the implementation of non-uniform grids for mag- 
netic field models, of UHE nuclei and secondary neutrinos and EM particles 
from their interactions, of inhomogeneous low energy target photon back- 
grounds for the UHE nuclei and EM cascade interactions, and of hadronic 
interactions with the baryon gas in dense parts of the large scale structure. 
Finally, interactions of UHE neutrinos with relic neutrinos of arbitrary mass 
and clustering properties could also be implemented, including the resulting 
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secondary particles. 
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